/** File Description
 *
 * Windbull (Wind statistical analyser)
 * Copyright (C) 2009 Ismail SEZEN
 * -------------------------------------------------------------------------------------------
 * This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 * -------------------------------------------------------------------------------------------
 * Project Name: Windbull
 * Files(4) : Windbull.c, Windbull.h, WindbullParser.c, WindbullParser.h
 * AUTHOR: Ismail SEZEN
 * e-mail: sezenismail@gmail.com
 * 13-05-2009 [12:20]
 * -------------------------------------------------------------------------------------------
 * Created by Code:Blocks 8.02
 */

/**Include Directives */
//{
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "WindbullParser.h"
//}

/** Defines */
//{
#define HOUR 24
#define RAD(val) ((M_PI/180) * (val)) //Converts degree to Radian
#define DEG(val) ((180/M_PI) * (val)) //Converts the radian to degree
#define RX(s,d) ((s) * cos(RAD(360.0-(((d)-5)*22.5)))) //Returns the x-component of wind vector
#define RY(s,d) ((s) * sin(RAD(360.0-(((d)-5)*22.5)))) //Returns the y-component of wind vector
#define ROUND(v) ( (v>=0) ? (int)(v+0.5) : (int)(v-0.5) ) //Returns the rounded value of double precision value
#define AIRDENSITY(h) ( (1.225 - (1.194e-04) * (h)) ) //Returns the density of th air at specified height
//}


/** Private*/
//{

/** Variables */
//{
int _wb_is_initialized = 0; // false
//}

/** Routines */

//{ Basic Routines

int compare_int(const void *i,const void *j){
	return *(int *)i - *(int *)j;
}

int get_class_interval(int data_count){
	if(data_count <= 0) return -1;
	double m = 1 + 3.3 * log10((double)data_count);
	int result = ROUND(m);
	return result;
}

int get_date_interval(char *pData,int buffer_size, int length,int *year,int *month,int *start,int *end){
	int id,y,m,d,si;
	char *format_temp = "%d;%d;%d;%d";
	(*start) = (*end) = 0;
	register int i;
	for(i=0;i<length;i++){
		char *strLine = (pData + (i * buffer_size));
		//printf("%s\n",strLine);
		//if(i==2) exit(-1);
		if( (si = sscanf(strLine, format_temp, &id, &y,&m,&d)) == 4){//Date handled suffesfully
			if(i == 0) *start = y;
			if(year == NULL){
				if(*start > y) *start = y;
				if(*end < y) *end = y;
			}else{
				if(month == NULL){
					if(*year == y){
						if(*start > m) *start = m;
						if(*end < m) *end = m;
					}
				}else{
					if(*year == y && *month == m){
						if(*start > d) *start = d;
						if(*end < d) *end = d;
					}
				}
			}
		}else{
			return -1;
		}
	}
	return 0;
}

double wind_speed(float *value, float *height){
	if(*height < 0){
		printf("Height must be positive for wind speed average. Application terminated! %g\n",*height);
		exit(1);
	}
	return *value;
}

double wind_power(float *value, float *height){
	double result = -1;
	if(*height > 0){
		double density = AIRDENSITY(*height);
		result = ( 0.5 * density * pow((*value),3) );
	}else{
		printf("Height must be positive for wind power. Application terminated!\n");
		exit(1);
	}
	return result;
}

//}

float get_max_speed(float *p_speed,int length){
	register int i;
	float result = *(p_speed);
	for(i=1;i<length;i++)
		if(*(p_speed + i) > result) result = *(p_speed + i);
	return result;
}

float get_min_speed(float *p_speed,int length){
	register int i;
	int result = *(p_speed);
	for(i=0;i<length;i++)
		if(*(p_speed + i) < result) result = *(p_speed + i);
	return result;
}

float speed_median(float *p_speed,int length){
	float *p_speed_copy;
	if((p_speed_copy = malloc(sizeof(float) * length)) == NULL){
		printf("Allocation Error\n");
		exit(1);
	}
	//Copy real array to other
	register int i;
	for(i=0;i<length;i++) *(p_speed_copy + i) = *(p_speed + i);
	//
	//Quicksort algorithm for median calculation
	qsort(p_speed_copy,length,sizeof(int),compare_int);
	//
	float median = (length&1) ? *(p_speed_copy + ((length+1)/2)) : ( (*(p_speed_copy + (((length)/2)-1)) + *(p_speed_copy + (((length)/2)+1)) ) / 2 );
	free(p_speed_copy);
	return median;
}

float speed_mode(float *p_speed,int length){
	float max = get_max_speed(p_speed,length);
	int mode_length = max * 10;
	int *p_mode_array;
	if((p_mode_array = malloc(sizeof(int) * mode_length)) == NULL){
		printf("Allocation Error\n");
		exit(1);
	}
	//
	register int i;
	//Reset the mode array
	for(i=0;i<mode_length;i++) *(p_mode_array) = 0;
	//
	//count the same speed values
	for(i=0;i<length;i++){
		int index = *(p_speed + i) * 10;
		(*(p_mode_array + index))++;
	}
	//Find most frequent index
	int most_freq = *(p_mode_array);
	float most_freq_val = 0.f;
	for(i=1;i<max;i++){
		if(*(p_mode_array + i) > most_freq){
			most_freq = *(p_mode_array + i);
			most_freq_val = (float)i/10.f;
		}
	}
	//
	free(p_mode_array);
	return most_freq_val;
}

float speed_avg(float *p_speed,int length){
	register int i;
	float total = 0.0;
	for(i=0;i<length;i++)
		total += *(p_speed + i);
	return total / length;
}

void wind_avg_hourly(float *p_speed, int length, double *hours, float *height,  double (*operation)(float *, float *) ){
	register int i,j;
	int hours_count = 24;
	int line_count = length / hours_count;
	//Reset hours array
	for(i=0;i<hours_count;i++) *(hours + i) = 0.0;
	//Sum each hours value
	for(i=0;i<hours_count;i++){
		for(j=0;j<line_count;j++){
			float *value  = ((float *)p_speed +(j * hours_count) + i);
			*(hours + i) += ((*operation)(value, height));
		}
	}
	//Division
	for(i=0;i<hours_count;i++)
		*(hours + i) = *(hours + i) / (double)line_count;
}

double variance(float *p_speed, int length){
	float avg = speed_avg(p_speed,length);
	float total = 0;
	register int i;
	for(i=0;i<length;i++)
		total += pow((*(p_speed + i) - avg),2);
	return ( total / (double)(length -1) );
}

void variance_hourly(float *p_speed, int length, float *height, double *p_avg_hourly_variance) {
	//
	register int i,j;
	int hours_count = 24;
	int line_count = length / hours_count;
	//
	double avg_hourly_speed[24];
	wind_avg_hourly(p_speed,length,&avg_hourly_speed[0], height, wind_speed);
	//
	//Reset hours array
	for(i=0;i<hours_count;i++) *(p_avg_hourly_variance + i) = 0.0;
	//Sum each hours value
	for(i=0;i<hours_count;i++){
		double total_for_var = 0.0;
		for(j=0;j<line_count;j++){
			float value  = *((float *)p_speed +(j * hours_count) + i);
			total_for_var += pow((value - avg_hourly_speed[i]),2);
		}
		*(p_avg_hourly_variance + i) = ( total_for_var / (double)(line_count -1) );
	}
}

double standart_deviation(float *p_speed,int length){
	return sqrt(variance(p_speed,length));
}

void standart_deviation_hourly(float *p_speed, int length, float *height, double *p_sd_hourly){
	variance_hourly(p_speed, length, height, &p_sd_hourly[0]);
	register int i;
	for(i=0;i<24;i++) *(p_sd_hourly+i) = sqrt(*(p_sd_hourly+i));
}

double autocorelation(float *p_speed,int length,int offset){
	if(offset > length/2) return -1;//offset can not greater than half length
	//
	float avg = speed_avg(p_speed,length);
	double v = variance(p_speed,length);
	double total = 0.0;
	//
	register int i;
	for(i=0;i<=(length -offset);i++){
		total += (*(p_speed +i) - avg) * (*(p_speed + i + offset) - avg);
	}
	double r = total / ( (length - offset) * v );
	return r;
}

void Wind_avg_monthly(char *pData, int buffer_size, float *p_speed, int length, double *avg_monthly,float *height,
										double (*operation)(float *, float *)){
	int hours_count = 24;
	int day_index = 0;
	int months_count = 12;
	double monthly_total = 0.0;
	int data_count[12] = {0,0,0,0,0,0,0,0,0,0,0,0};
	register int i,j,k,l;
	//--------------------
	//Reset hours array
	for(i=0;i<months_count;i++) *(avg_monthly + i) = 0.0;
	//--------------------
	int year_start = 0,year_end = 0;
	get_date_interval(pData,buffer_size,length,NULL,NULL,&year_start,&year_end);
	int year_count = year_end - year_start + 1;
	for(i=0;i<=year_count;i++){
		int month_start  = 0, month_end = 0;
		int year = year_start + i;
		get_date_interval(pData,buffer_size,length,&year,NULL,&month_start,&month_end);
		//int month_interval = month_end - month_start + 1;
		for(j=month_start-1;j<month_end;j++){
			int day_start  = 0, day_end = 0;
			int month = j +1;
			get_date_interval(pData,buffer_size,length,&year,&month,&day_start,&day_end);
			int day_interval = day_end - day_start + 1;
			data_count[j] += day_interval * hours_count;
			monthly_total = 0.0;
			for(k=day_start-1;k<day_end;k++){
				for(l=0;l<hours_count;l++){
					float *value  = ((float *)p_speed + (day_index * hours_count) + l);
					//printf("%g ",*value);
					//float *value  = ((float *)p_speed + (i * year_count) + (j * month_interval) + (k * day_interval) + l);
					monthly_total += ((*operation)(value, height));
				}
				day_index++;
			}
			*(avg_monthly+j) += monthly_total;
		}
	}
	//printf("day index = %d\n",day_index);
	//--------------------
	for(i=0;i<months_count;i++){
		*(avg_monthly+i) = *(avg_monthly+i) / (double) data_count[i];
	}
}

void Wind_avg_seasonal(char *pData, int buffer_size, float *p_speed, int length, double *avg_seasonal,float *height,
										double (*operation)(float *, float *)){
	int hours_count = 24;
	int day_index = 0;
	int season_count = 4;
	int data_count[4] = {0,0,0,0};
	register int i,j,k,l;
	//--------------------
	//Reset hours array
	for(i=0;i<season_count;i++) *(avg_seasonal + i) = 0.0;
	//--------------------
	int year_start = 0,year_end = 0;
	get_date_interval(pData,buffer_size,length,NULL,NULL,&year_start,&year_end);
	int year_count = year_end - year_start + 1;
	for(i=0;i<=year_count;i++) {
		int month_start  = 0, month_end = 0;
		int year = year_start + i;
		get_date_interval(pData,buffer_size,length,&year,NULL,&month_start,&month_end);
		//int month_interval = month_end - month_start + 1;
		for(j=month_start-1;j<month_end;j++) {
			int day_start  = 0, day_end = 0;
			int month = j +1;
			get_date_interval(pData,buffer_size,length,&year,&month,&day_start,&day_end);
			int day_interval = day_end - day_start + 1;
			double monthly_data_count = day_interval * hours_count;
			double monthly_total = 0.0;
			for(k=day_start-1;k<day_end;k++) {
				for(l=0;l<hours_count;l++) {
					float *value  = ((float *)p_speed + (day_index * hours_count) + l);
					//float *value  = ((float *)p_speed + (i * year_count) + (j * month_interval) + (k * day_interval) + l);
					monthly_total += ((*operation)(value, height));
				}
				day_index++;
			}
			//
			switch (j) {
				case 0:case 1:case 11:
					//printf("Winter\t");
					*(avg_seasonal) += monthly_total;
					data_count[0] += monthly_data_count;
					break;
				case 2:case 3:case 4:
					//printf("Spring\t");
					*(avg_seasonal + 1) += monthly_total;
					data_count[1] += monthly_data_count;
					break;
				case 5:case 6:case 7:
					//printf("Summer\t");
					*(avg_seasonal + 2) += monthly_total;
					data_count[2] += monthly_data_count;
					break;
				case 8:case 9:case 10:
					//printf("Autumn\t");
					*(avg_seasonal + 3) += monthly_total;
					data_count[3] += monthly_data_count;
					break;
			}
			//
		}
	}
	//--------------------
	for(i=0;i<season_count;i++){
		*(avg_seasonal+i) = *(avg_seasonal+i) / (double) data_count[i];
	}
}

void Wind_avg_monthly_day_night(char *pData, int buffer_size, float *p_speed, int length, double *avg_monthly_day, double *avg_monthly_night, float *height,
										double (*operation)(float *, float *)){
	int hours_count = 24;
	int day_index = 0;
	int months_count = 12;
	int data_count[12] = {0,0,0,0,0,0,0,0,0,0,0,0};
	register int i,j,k,l;
	//--------------------
	//Reset hours array
	for(i=0;i<months_count;i++){
		 *(avg_monthly_day + i) = 0.0;
		 *(avg_monthly_night + i) = 0.0;
	}
	//--------------------
	int year_start = 0,year_end = 0;
	get_date_interval(pData,buffer_size,length,NULL,NULL,&year_start,&year_end);
	int year_count = year_end - year_start + 1;
	for(i=0;i<=year_count;i++) {
		int month_start  = 0, month_end = 0;
		int year = year_start + i;
		get_date_interval(pData,buffer_size,length,&year,NULL,&month_start,&month_end);
		//int month_interval = month_end - month_start + 1;
		for(j=month_start-1;j<month_end;j++) {
			int day_start  = 0, day_end = 0;
			int month = j +1;
			get_date_interval(pData,buffer_size,length,&year,&month,&day_start,&day_end);
			//int day_interval = day_end - day_start + 1;
			double monthly_total_day = 0.0;
			double monthly_total_night = 0.0;
			for(k=day_start-1;k<day_end;k++) {
				for(l=0;l<hours_count;l++) {
					float *value  = ((float *)p_speed + (day_index * hours_count) + l);
					if( j > 8 || j < 3 ){ //Winter
							if( l > 6 || l < 18 ){ // 7 - 17
								data_count[j]++;
								monthly_total_day += ((*operation)(value, height));
							}
							if( l < 7 || l > 17 ){ // 18 - 06
								data_count[j]++;
								monthly_total_night += ((*operation)(value, height));
							}
					}else{ // Summer
							if( l > 5 || l < 19 ){ // 06 - 18
								data_count[j]++;
								monthly_total_day += ((*operation)(value, height));
							}
							if( l > 18 || l < 6 ){ //19 - 05
								data_count[j]++;
								monthly_total_night += ((*operation)(value, height));
							}
					}
				}
			}
			day_index++;
			*(avg_monthly_day + j) += monthly_total_day;
			*(avg_monthly_night + j) += monthly_total_night;
		}
	}
	//--------------------
	for(i = 0; i < months_count; i++){
		*(avg_monthly_day + i) = *(avg_monthly_day + i) / (double) data_count[i];
		*(avg_monthly_night + i) = *(avg_monthly_night + i) / (double) data_count[i];
	}
}

void variance_monthly(char *pData, int buffer_size, float *p_speed, int length, float *height,  double *p_array_12){
	double avg_monthly[12];
	Wind_avg_monthly(pData, buffer_size, p_speed, length,&avg_monthly[0],height,wind_speed);
	int hours_count = 24;
	int day_index = 0;
	int months_count = 12;
	double monthly_total[12] = {0,0,0,0,0,0,0,0,0,0,0,0};
	int data_count[12] = {0,0,0,0,0,0,0,0,0,0,0,0};
	register int i,j,k,l;
	//--------------------
	//Reset hours array
	for(i=0;i<months_count;i++) *(p_array_12 + i) = 0.0;
	//--------------------
	int year_start = 0,year_end = 0;
	get_date_interval(pData,buffer_size,length,NULL,NULL,&year_start,&year_end);
	int year_count = year_end - year_start + 1;
	//printf("year count : %d\n",year_count);
	for(i=0;i<year_count;i++){
		int month_start  = 0, month_end = 0;
		int year = year_start + i;
		get_date_interval(pData,buffer_size,length,&year,NULL,&month_start,&month_end);
		//printf("month_count = %d\n",(month_end - month_start + 1));
		//int month_interval = month_end - month_start + 1;
		for(j=month_start-1;j<month_end;j++){
			int day_start  = 0, day_end = 0;
			int month = j +1;
			get_date_interval(pData,buffer_size,length,&year,&month,&day_start,&day_end);
			int day_interval = day_end - day_start + 1;
			data_count[j] += day_interval * hours_count;
			//monthly_total = 0.0;
			//printf("%d %d\n",day_start,day_end);
			for(k=day_start-1;k<day_end;k++){
				for(l=0;l<hours_count;l++){
					float *value  = ((float *)p_speed + (day_index * hours_count) + l);
					monthly_total[j] += pow(((*value) - avg_monthly[j]),2);
				}
				day_index++;
			}
		}
	}
	/*--------------------*/
	for(i=0;i<months_count;i++){
		*(p_array_12+i) = monthly_total[i] / ((double) (data_count[i]-1));
	}
}

void sd_monthly(char *pData, int buffer_size, float *p_speed, int length, float *height,  double *p_array_12){
	double var_array[12];
	/* Get the variance values of the each months*/
	variance_monthly(pData, buffer_size, p_speed, length, height, &var_array[0]);
	register int i;
	/* Calculate the square root of each variance belong to the month*/
	for(i=0;i<12;i++) *(p_array_12+i) = sqrt(var_array[i]);
}

///----------------------------------------------------------------------------------------------------------------------------------------
double power_coefficient_length(float z0){
	return (0.096f * log10(z0) + 0.016f * pow(log10(z0),2) + 0.24f);
}

double power_coefficient_speed_and_height(float Uref, float Zref){
	if(Uref == 0) Uref = 0.0001f;
	return ( (0.37f - 0.088 * log (Uref)) / (1.0 - 0.088f * log(Zref/10) ) );
}

double power_coefficient_speed_and_length(float Uref, float Z0){
	if(Uref == 0) Uref = 0.0001f;
	return Z0 * (1 - 0.55f * log(Uref));
}
///----------------------------------------------------------------------------------------------------------------------------------------

double extraplation_height(float Uref, float Zref, float Z,double a){
	return Uref * pow((Z / Zref),a);
}

void extrapolate_speed(float *p_speed,float *p_speed_extrapolated, int length, float Zref, float Z, float z0){
	//float *p_speed_copy;
	printf("Zref = %g Z = %g \n",Zref,Z);
	if(p_speed_extrapolated == NULL){
		if((p_speed_extrapolated = (float *) malloc(sizeof(float) * length)) == NULL){
			printf("Allocation Error\n");
			exit(1);
		}
	}
	/*Copy real array to other*/
	register int i;
	for(i=0;i<length;i++){
		float Uref = *(p_speed+i);
		double a = power_coefficient_speed_and_height(Uref, Zref);
		//double a = power_coefficient_speed_and_length(Uref, z0);
		//double a = power_coefficient_length(z0);
		*(p_speed_extrapolated + i) = extraplation_height(Uref, Zref, Z, a);
		//if((i % 24) == 0) printf("\nL(%d) ",(i/24));
		//printf("%g ",*(p_speed_extrapolate + i));
		//printf("%d =%g\n",i,*(p_speed_extrapolate + i));
	}
}

double energy_pattern_factor(float *p_speed,int length){
	int total3 = 0;
	register int i;
	for(i=0;i<length;i++){
		total3 += pow(*(p_speed+i),3);
	}
	double avg3 = total3 / length;
	float avg = speed_avg(p_speed,length);
	return avg3 / pow(avg,3);
}

void speed_frequency(float *p_speed,int length, double *p_ui, int *p_mi, double *p_pui, double *p_Fui, float *class_bin, int *class_interval){
	int w = get_class_interval(length);
	*class_interval = w;
	float max = get_max_speed(p_speed,length);
	float bin = max  / (float) (w);
	*class_bin = bin;
	//Alloc memory for arrays
	//double *p_ui;
	if((p_ui = malloc(sizeof(double) * w)) == NULL){
		printf("Allocation Error\n");exit(1);
	}
	//int *p_mi;
	if((p_mi = malloc(sizeof(int) * w)) == NULL){
		printf("Allocation Error\n");exit(1);
	}
	//double *p_pui;
	if((p_pui = malloc(sizeof(double) * w)) == NULL){
		printf("Allocation Error\n");exit(1);
	}
	//double *p_Fui;
	if((p_Fui = malloc(sizeof(double) * w)) == NULL){
		printf("Allocation Error\n");exit(1);
	}
	//
	register int i,j;
	//Prepare p_mi array for first use
	for(i=0;i<w;i++) *(p_mi+i) = 0;
	//Calculate ui values
	for(i=0;i<w;i++)
		*(p_ui + i) = (i * bin) + (bin/2);
	//Calculate mi values
	for(i=0;i<length;i++){
		int counted = 1;
		for(j=0;j<w;j++){
			if( ( *(p_speed+i) >= (j * bin) ) && ( *(p_speed +i) < ( (j+1)*bin ) ) ) {
				(*(p_mi + j))++;
				counted = 0;
			}
		}
		if(counted == 1){
			if(*(p_speed + i) >= max){
				(*(p_mi + w -1))++;
			}
		}
	}
	//Calculate pui probability values
	for(i=0;i<w;i++){
		*(p_pui + i) = *(p_mi+i) / (float)length;
	}
	//
	//Calculate Fui Cumulative Probability values
	double total = 0.0;
	for(i=0;i<w;i++){
		total += *(p_mi + i);
		*(p_Fui + i) = ( total )/(float)length;
	}
	//Print results
	printf("\nBin value = %g\n",bin);
	printf("Class count = %d\n",w);
	printf("i\tinterval\t\t\tui\t\tmi\tp(ui)\t\tF(ui)\n");
	printf("----------------------------------------------------------------------------------------------\n");
	for(i=0;i<w;i++){
		printf("%d\t%e-%e\t%e\t%d\t%e\t%e\n",i,(i * bin),((i+1)*bin),(*(p_ui+i)),*(p_mi+i),*(p_pui+i),*(p_Fui+i));
	}
//	//Free memory
//	free(p_ui);
//	free(p_Fui);
//	free(p_mi);
//	free(p_pui);
}

inline int Teta(double value){
	return (value>=0)?1:0;
}

inline double Ti(float *p_array_speed,double ixdelta_V, int N){
	register int j;
	int teta_total = 0;
	for(j=0;j<N;j++){
		teta_total+=Teta(*(p_array_speed+j) - ixdelta_V);
	}
	return teta_total / (double)N;
}

int SDC(float *p_array_speed, int length, double** p_array_Ti, double** p_array_Wt){
	int D = get_class_interval(length);
	float max = get_max_speed(p_array_speed,length);
	float delta_V = max  / (float) (D);
	D++;
	double* p_Ti = NULL;
	double *p_Wt = NULL;
	if((p_Ti = (double *) malloc(sizeof(double) * D)) == NULL){
		printf("Allocation Error\n");exit(1);
	}
	if((p_Wt = (double *) malloc(sizeof(double) * D)) == NULL){
		printf("Allocation Error\n");exit(1);
	}
	register int i;
	for(i=0;i<D;i++){
		double ixdelta_V = (i * delta_V);
		double Ti_result = Ti(p_array_speed,ixdelta_V, length)*100;
		*(p_Wt+i) = ixdelta_V;
		*(p_Ti+i) =Ti_result;
	}
	*p_array_Ti = p_Ti;
	*p_array_Wt = p_Wt;
	return D;
}

void wind_direction_frequencies(int *p_direction,int data_count, int *p_frequency){
	register int i;
	for(i=0;i<=16;i++) *(p_frequency+i) = 0; //Reset the array
	for(i=0;i<data_count;i++){
		int direc = *(p_direction + i);
		(*(p_frequency + direc))++;
	}
}

void wind_speed_frequencies(int *p_direction, float *p_speed, int length,double *p_frequency){
	register int i;
	int data_count[17];
	/*Reset the array*/
	for(i=0;i<=16;i++){
		 *(p_frequency+i) = 0;
		 data_count[i] = 0;
	}
	/**/
	for(i=0;i<length;i++){
		float val = *(p_speed+i);
		int direc = *(p_direction + i);
		(*(p_frequency + direc))+=val;
		data_count[direc]++;
	}
	/**/
	for(i=0;i<=16;i++){
		*(p_frequency + i) = *(p_frequency + i) / (double)data_count[i];
	}
}

void wind_direction_frequencies_print(int *p_direction,int data_count){
	register int i;
	int dir_freq[16] ={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
	wind_direction_frequencies(p_direction,data_count,&dir_freq[0]);
	printf("Directions:\t");
	for(i=0;i<16;i++) printf("%d\t",(i+1));
	printf("\n");
	printf("----------\t");
	for(i=0;i<16;i++) printf("------\t");
	printf("\nCounts:\t\t");
	int total = 0;
	for(i=0;i<16;i++){
		total+=dir_freq[i];
		printf("%d\t",dir_freq[i]);
	}
	printf("\n");
	printf("Total defined directions count = %d\n",total);
	printf("Total undefined directions count = %d\n",(data_count - total) );
}

double weibull_k_param(double sd, float avg_speed){
	return pow((sd/(double)avg_speed),(-1.086));
}

double gamma(double k){
	double n = 1/k;
	return pow(n,n) * exp(-n) * sqrt(2 * M_PI * n) * ( 1 + (1 / (12 * n)) + (1 / (288 * n * n)) - (139.0/(51840 * n *n * n)) );
}

inline void weibull_k_param_by_sd_and_avg(double *p_array_sd, double *p_array_avg, int length, double *p_array_result){
	register int i;
	for(i=0;i<length;i++){
		*(p_array_result+i) = pow( (p_array_sd[i] / p_array_avg[i]),(-1.086) );
	}
}

void weibull_k_param_monthly(float *pData, int buffer_size, float *p_speed, int length, float *height, double *k_array_12){
	/*Calculate the standart deviation for each month*/
	double sd_array_12[12];
	sd_monthly(pData, buffer_size, p_speed, length, height, &sd_array_12[0]);
	/*Calculate the avg speed for each month*/
	double avg_monthly[12];
	Wind_avg_monthly(pData, buffer_size, p_speed, length,&avg_monthly[0],height,wind_speed);
	/*Calculate the k parameter for each month*/
	weibull_k_param_by_sd_and_avg(&sd_array_12[0], &avg_monthly[0], 12, k_array_12);
}

double weibull_c_param(float avg_speed, double k){
return avg_speed / gamma(k);
}

void weibull_c_param_monthly(double *p_array_avg_speed, double *p_array_k_param, int length, double *p_array_result){
	register int i;
	for(i=0;i<length;i++){
		p_array_result[i] = ((double)p_array_avg_speed[i]) / gamma(p_array_k_param[i]);
	}
}

///----------------------------------------------------------------------------------------------------------------------------------------
///Lambert Persistance  Routines
///----------------------------------------------------------------------------------------------------------------------------------------

double hourly(double *p_array,int j_length, int i_hours_count){
	int hours_count = 24;
	double *value  = ((double *)p_array +(j_length * hours_count) + i_hours_count);
	return *value;
}

void get_lambert_persistance_value_detailed(float *p_speed,int *p_direction,int length,
																		float *total_wind_speed,
																		float *avg_wind_speed, float *total_Rx,
																		float *total_Ry,float *Rx,float *Ry,
																		float *R,float *p,float *B){
		register int i;
		//Vectoral Summation of direction and wind speeds and
		//total wind summation for each wind value
		for(i = 0; i < length; i++) {
			float speed_val  = *(p_speed + i);
			float direction_val  = *(p_direction + i);
			//
			(*total_wind_speed) += speed_val;
			(*total_Rx) += RX(speed_val, direction_val);
			(*total_Ry) += RY(speed_val, direction_val);
		}
		//Lambert interval and latest Results
		(*avg_wind_speed) = (*total_wind_speed) / (float)(length);
		(*Rx) = ( (*total_Rx) * (-1)) / (float)(length);
		(*Ry) = ( (*total_Ry) * (-1)) / (float)(length);
		(*R) = sqrt((pow( (*Rx), 2) + pow( (*Ry), 2)));
		(*p) = (*R) / (*avg_wind_speed);
		(*B) = atan((double)((*Ry)/(*Rx)));
}

void get_lambert_persistance_value_detailed2(float *p_speed,int *p_direction,int length,
																		float *total_wind_speed,
																		float *avg_wind_speed, float *total_Rx,
																		float *total_Ry,float *Rx,float *Ry,
																		float *R,float *p,float *B){
		register int i;
		float alfa = RAD(22.5f);
		float Rx_total_direc_1 = 0.0f, Rx_total_direc_2 = 0.0f, Rx_total_direc_3 = 0.0f;
		float Ry_total_direc_1 = 0.0f, Ry_total_direc_2 = 0.0f, Ry_total_direc_3 = 0.0f;
		//Vectoral Summation of direction and wind speeds and
		//total wind summation for each wind value
		for(i = 0; i < length; i++) {
			float speed_val  = *(p_speed + i);
			int direction_val  = *(p_direction + i);
			//
			(*total_wind_speed) += speed_val;
			//For Rx
			switch(direction_val){
				case 13: // W
					(*total_Rx) += speed_val;break;
				case 5: // E
					(*total_Rx) -= speed_val;break;
				case 10:case 16: // SSW and NNW
					Rx_total_direc_1 += speed_val;break;
				case 8:case 2: // SSE and NNE
					Rx_total_direc_1 -= speed_val;break;
				case 11:case 15: // SW and NW
					Rx_total_direc_2 += speed_val;break;
				case 3:case 7: // NE and SE
					Rx_total_direc_2 -= speed_val;break;
				case 12:case 14: // WSW and WNW
					Rx_total_direc_3 += speed_val;break;
				case 4:case 6: //ENE and ESE
					Rx_total_direc_3 -= speed_val;break;
			}
			//For Ry
			switch(direction_val){
				case 9: // S
					(*total_Ry) += speed_val;break;
				case 1: // N
					(*total_Ry) -= speed_val;break;
				case 8:case 10: // SSE and SSW
					Ry_total_direc_1 += speed_val;break;
				case 2:case 16: // NNE NNW
					Ry_total_direc_1 -= speed_val;break;
				case 7:case 11: // SE and SW
					Ry_total_direc_2 += speed_val;break;
				case 3:case 15: // NE and NW
					Ry_total_direc_2 -= speed_val;break;
				case 12:case 6: // WSW and ESE
					Ry_total_direc_3 += speed_val;break;
				case 14:case 4: // WNW and ENE
					Ry_total_direc_3 -= speed_val;break;
			}
		}
		(*total_Rx) += (Rx_total_direc_1 * sin(alfa)) + ( Rx_total_direc_2 * sin(2 * alfa) ) + ( Rx_total_direc_3 * sin(3 * alfa) );
		(*total_Ry) += (Ry_total_direc_1 * cos(alfa)) + ( Ry_total_direc_2 * cos(2 * alfa) ) + ( Ry_total_direc_3 * cos(3 * alfa) );
		//Lambert interval and latest Results
		(*avg_wind_speed) = (*total_wind_speed) / (float)(length);
		(*Rx) = ( (*total_Rx)) / (float)(length);
		(*Ry) = ( (*total_Ry)) / (float)(length);
		(*R) = sqrt((pow( (*Rx), 2) + pow( (*Ry), 2)));
		(*p) = (*R) / (*avg_wind_speed);
		(*B) = atan((double)((*Rx)/(*Ry)));
}

double lambert_persistance_value(float *p_speed,int *p_direction,int length){
	float total_Rx=0.0,total_Ry=0.0,total_wind_speed=0.0,avg_wind_speed =0.0;
	float Rx, Ry, R, p,B;
	get_lambert_persistance_value_detailed(p_speed,p_direction, length,
																&total_wind_speed,&avg_wind_speed,
																&total_Rx,&total_Ry,&Rx,&Ry,&R,&p,&B);
	return p;
}

void lambert_persistant_value_hourly(float *p_speed,int *p_direction,int length, float *p_persistant_values){

	register int i,j;
	int hours_count = 24;
	int line_count = length / hours_count;
	//Reset hours array
	for(i = 0; i < hours_count; i++) *(p_persistant_values + i) = 0.0f;
	//
	//Calculate each hour's value
	for(i = 0; i < hours_count; i++) {
		float total_Rx=0.0, total_Ry=0.0, total_wind_speed=0.0;
		float Rx_total_direc_1 = 0.0f, Rx_total_direc_2 = 0.0f, Rx_total_direc_3 = 0.0f;
		float Ry_total_direc_1 = 0.0f, Ry_total_direc_2 = 0.0f, Ry_total_direc_3 = 0.0f;
		for(j = 0; j < line_count; j++) {
			float speed_val  = *((float *)p_speed +(j * hours_count) + i);
			int direction_val  = *((int *)p_direction +(j * hours_count) + i);
			total_wind_speed += speed_val;
			total_Rx += RX(speed_val, direction_val);
			total_Ry += RY(speed_val, direction_val);
		}
		//Lambert interval and latest Results
		float avg_wind_speed = total_wind_speed / (float)(line_count);
		float Rx = total_Rx / (float)(line_count);
		float Ry = total_Ry / (float)(line_count);
		float R = sqrt((pow(Rx,2) + pow(Ry,2)));
		*(p_persistant_values+i) = R / avg_wind_speed;
	}
}

void lambert_persistant_value_monthly(char *pData, int buffer_size, float *p_speed,int *p_direction,int length, float *p_persistant_values){
	float alfa = RAD(22.5f);
	int hours_count = 24;
	int months_count = 12;
	int data_count[12];
	float total_Rx[12], total_Ry[12], total_wind_speed[12];
//	float Rx_total_direc_1[12], Rx_total_direc_2[12], Rx_total_direc_3[12];
//	float Ry_total_direc_1[12], Ry_total_direc_2[12], Ry_total_direc_3[12];
	int day_index = 0;
	register int i,j,k,l;
	//--------------------
	//Reset and init arrays
	for(i=0;i<months_count;i++) {
		*(p_persistant_values + i) = 0.0f;
		data_count[i] = 0;
		total_Rx[i] = 0.0f;
		total_Ry[i] = 0.0f;
		total_wind_speed[i] = 0.0f;
//		Rx_total_direc_1[i] = 0.0f;
//		Rx_total_direc_2[i] = 0.0f;
//		Rx_total_direc_3[i] = 0.0f;
//		Ry_total_direc_1[i] = 0.0f;
//		Ry_total_direc_2[i] = 0.0f;
//		Ry_total_direc_3[i] = 0.0f;
	}
	//--------------------
	int year_start = 0,year_end = 0;
	get_date_interval(pData,buffer_size,length,NULL,NULL,&year_start,&year_end);
	int year_count = year_end - year_start + 1;
	for(i = 0; i <= year_count; i++){
		int month_start  = 0, month_end = 0;
		int year = year_start + i;
		get_date_interval(pData,buffer_size,length,&year,NULL,&month_start,&month_end);
		int month_interval = month_end - month_start + 1;
		for(j = month_start - 1; j < month_end; j++) {
			int day_start  = 0, day_end = 0;
			int month = j +1;
			get_date_interval(pData,buffer_size,length,&year,&month,&day_start,&day_end);
			int day_interval = day_end - day_start + 1;
			data_count[j] += day_interval * hours_count;
			//printf("%d Data count for month = %d %d\n",(j+1),day_interval,data_count[j]);
			for(k = day_start - 1; k < day_end; k++){
				for(l = 0; l < hours_count; l++) {
					//float speed_val  = *((float *)p_speed + (i * year_count) + (j * month_interval) + (k * day_interval) + l);
					//loat direction_val  = *((float *)p_direction + (i * year_count) + (j * month_interval) + (k * day_interval) + l);
					float speed_val  = *((float *)p_speed + (day_index * hours_count) + l);
					int direction_val  = *((int *)p_direction + (day_index * hours_count) + l);
					//
					total_wind_speed[j] += speed_val;
					total_Rx[j] += RX(speed_val, direction_val);
					total_Ry[j] += RY(speed_val, direction_val);
//					//For Rx
//					switch(direction_val){
//						case 13: // W
//							total_Rx[j] += speed_val;break;
//						case 5: // E
//							total_Rx[j] -= speed_val;break;
//						case 10:case 16: // SSW and NNW
//							Rx_total_direc_1[j] += speed_val;break;
//						case 8:case 2: // SSE and NNE
//							Rx_total_direc_1[j] -= speed_val;break;
//						case 11:case 15: // SW and NW
//							Rx_total_direc_2[j] += speed_val;break;
//						case 3:case 7: // NE and SE
//							Rx_total_direc_2[j] -= speed_val;break;
//						case 12:case 14: // WSW and WNW
//							Rx_total_direc_3[j] += speed_val;break;
//						case 4:case 6: //ENE and ESE
//							Rx_total_direc_3[j] -= speed_val;break;
//					}
//					//For Ry
//					switch(direction_val){
//						case 9: // N
//							total_Ry[j] += speed_val;break;
//						case 1: // N
//							total_Ry[j] -= speed_val;break;
//						case 8:case 10: // SSE and SSW
//							Ry_total_direc_1[j] += speed_val;break;
//						case 2:case 16: // NNE NNW
//							Ry_total_direc_1[j] -= speed_val;break;
//						case 7:case 11: // SE and SW
//							Ry_total_direc_2[j] += speed_val;break;
//						case 3:case 15: // NE and NW
//							Ry_total_direc_2[j] -= speed_val;break;
//						case 12:case 6: // WSW and ESE
//							Ry_total_direc_3[j] += speed_val;break;
//						case 14:case 4: // WNW and ENE
//							Ry_total_direc_3[j] -= speed_val;break;
//					}
				}
				day_index++;
			}
		}
	}
	//--------------------
	for(i = 0; i < months_count; i++){
//		total_Rx[i] += (Rx_total_direc_1[i] * sin(alfa)) + ( Rx_total_direc_2[i] * sin(2 * alfa) ) + ( Rx_total_direc_3[i] * sin(3 * alfa) );
//		total_Ry[i] += (Ry_total_direc_1[i] * cos(alfa)) + ( Ry_total_direc_2[i] * cos(2 * alfa) ) + ( Ry_total_direc_3[i] * cos(3 * alfa) );
		//Lambert interval and latest Results
		float avg_wind_speed = total_wind_speed[i] / ((float) data_count[i]);
		float Rx = (total_Rx[i] * (-1)) / ((float) data_count[i]);
		float Ry = (total_Ry[i] * (-1)) / ((float) data_count[i]);
		float R = sqrt((pow(Rx,2) + pow(Ry,2)));
		*(p_persistant_values+i) = R / avg_wind_speed;
	}
}

void lambert_persistant_value_seasonal(char *pData, int buffer_size, float *p_speed,int *p_direction,int length, float *p_persistant_values){
	int hours_count = 24;
	int day_index = 0;
	int season_count = 4;
	int data_count[4];
	//
	float total_Rx[4], total_Ry[4], total_wind_speed[4];
	register int i,j,k,l;
	//--------------------
	//Reset and init arrays
	for(i=0;i<season_count;i++) {
		*(p_persistant_values + i) = 0.0f;
		data_count[i] = 0.0;
		total_Rx[i] = 0.0f;
		total_Ry[i] = 0.0f;
		total_wind_speed[i] = 0.0f;
	}
	//--------------------
	int year_start = 0,year_end = 0;
	get_date_interval(pData,buffer_size,length,NULL,NULL,&year_start,&year_end);
	int year_count = year_end - year_start + 1;
	for(i=0;i<=year_count;i++){
		int month_start  = 0, month_end = 0;
		int year = year_start + i;
		get_date_interval(pData,buffer_size,length,&year,NULL,&month_start,&month_end);
		int month_interval = month_end - month_start + 1;
		for(j=month_start-1;j<month_end;j++){
			int day_start  = 0, day_end = 0;
			int month = j +1;
			get_date_interval(pData,buffer_size,length,&year,&month,&day_start,&day_end);
			int day_interval = day_end - day_start + 1;
			double monthly_data_count = day_interval * hours_count;
			double monthly_total_wind_speed  = 0.0;
			double monthly_Rx = 0.0;
			double monthly_Ry = 0.0;
			for(k=day_start-1;k<day_end;k++) {
				for(l=0;l<hours_count;l++) {
					float speed_val  = *((float *)p_speed + (day_index * hours_count) + l);
					int direction_val  = *((int *)p_direction + (day_index * hours_count) + l);
					monthly_total_wind_speed += speed_val;
					monthly_Rx += RX(speed_val, direction_val);
					monthly_Ry += RY(speed_val, direction_val);
				}
				day_index++;
			}
			//
			switch (j) {
				case 0:case 1:case 11:
					//printf("Winter\t");
					total_wind_speed[0] += monthly_total_wind_speed;
					total_Rx[0] += monthly_Rx;
					total_Ry[0] += monthly_Ry;
					data_count[0] += monthly_data_count;
					break;
				case 2:case 3:case 4:
					//printf("Spring\t");
					total_wind_speed[1] += monthly_total_wind_speed;
					total_Rx[1] += monthly_Rx;
					total_Ry[1] += monthly_Ry;
					data_count[1] += monthly_data_count;
					break;
				case 5:case 6:case 7:
					//printf("Summer\t");
					total_wind_speed[2] += monthly_total_wind_speed;
					total_Rx[2] += monthly_Rx;
					total_Ry[2] += monthly_Ry;
					data_count[2] += monthly_data_count;
					break;
				case 8:case 9:case 10:
					//printf("Autumn\t");
					total_wind_speed[3] += monthly_total_wind_speed;
					total_Rx[3] += monthly_Rx;
					total_Ry[3] += monthly_Ry;
					data_count[3] += monthly_data_count;
					break;
			}
			//
		}
	}
	//--------------------
	for(i = 0; i < season_count; i++){
		//Lambert interval and latest Results
		float avg_wind_speed = total_wind_speed[i]/((double) data_count[i]);
		float Rx = (total_Rx[i] * (-1)) / ((double) data_count[i]);
		float Ry = (total_Ry[i] * (-1)) / ((double) data_count[i]);
		float R= sqrt((pow(Rx,2) + pow(Ry,2)));
		*(p_persistant_values+i) = R / avg_wind_speed;
		//*(p_persistant_values+i) = (R[i]);
	}
}

double lambert_persistance_value_print(float *p_speed,int *p_direction,int length){
		float total_Rx=0.0,total_Ry=0.0,total_wind_speed=0.0,avg_wind_speed;
		float Rx, Ry, R, p,B;
		get_lambert_persistance_value_detailed2(p_speed,p_direction, length,
																	&total_wind_speed,&avg_wind_speed,
																	&total_Rx,&total_Ry,&Rx,&Ry,&R,&p,&B);
		//
		printf("Lamber Persistant Method : \n");
		printf("----------------------------------\n");
		printf("Total Data Count = %d\n",length);
		printf("Total Wind Speed = %g\n",total_wind_speed);
		printf("Avg.Wind Speed = %g\n",avg_wind_speed);
		printf("\nRx Total = %e\tRy Total = %e\n",total_Rx,total_Ry);
		printf("Rx = %g Ry = %g\n",Rx,Ry);
		printf("R = %g\n",R);
		printf("Persistant value = %g\n",p);
		printf("Angle = %g\n",DEG(B));
		printf("----------------------------------\n");

		return p;
}

///----------------------------------------------------------------------------------------------------------------------------------------

//}

/** Public*/
//{

/** Enums */
enum WB_Directions_enum{CALM,N,NNE,NE,ENE,E,ESE,SE,SSE,S,SSW,SW,WSW,W,WNW,NW,NNW};
enum WB_Directions_enum WB_direction;
char  WB_Direction_to_string[17][5] = {"CALM","N","NNE","NE","ENE","E","ESE","SE","SSE","S","SSW","SW","WSW","W","WNW","NW","NNW"};
/** Variables */
//{
int WB_Total_Data_Count = 0;
float *WB_arr_speed = NULL;
float *WB_arr_speed_ext = NULL;
int *WB_arr_direction = NULL;
float WB_height = 10.0f;
float WB_height_old = 10.0f;
float WB_default_height = 10.0f;
float WB_z0 = 0.0001f; /*roughness length for flat surface*/
//
int WB_histo_class_count = 0;
int WB_histo_bin_value = 0;
double *WB_histo_ui = NULL;
int *WB_histo_mi = NULL;
double *WB_histo_pui = NULL;
double *WB_histo_fui = NULL;
//
float WB_speed_max_var = -1;
float WB_speed_min_var = -1;
float WB_speed_avg_var = -1;
float WB_speed_median_var = -1;
float WB_speed_mode_var = -1;
double WB_speed_sd_var = -1;
double WB_speed_variance_var = -1;
double WB_lambert_persistance_var = -1;
double WB_autocorrelation_var = -99;
double WB_energy_pattern_factor_var = -1;
double WB_weibull_k_param_var = -1;
double WB_weibull_c_param_var = -1;
//}

/** Routines */
//{
/**
 * Function: WB_Initialize
 * ----------------------------------
 * Initializes the Windbull. Reads and parses the file and creates needed arrays for process.
 * @param *fileName is the file will be read
 * returns : Zero (0) if  the operation is succeed otherwise returns subzero (-1).
*/
int WB_Initialize(char *fileName){
	printf("WB_Initialize (%s)\n",fileName);
	if(WBP_read_file(fileName) == 0){
		WB_Total_Data_Count = WBP_Length * HOUR;
		//Alloc enough memory for speed array
		if((WB_arr_speed = (float *) malloc(sizeof(float) * WB_Total_Data_Count)) == NULL){
			printf("Allocation Error for speed array\n");exit(1);
		}
		//Alloc enough space for direction array
		if((WB_arr_direction = (int *) malloc(sizeof(int) * WB_Total_Data_Count)) == NULL){
			printf("Allocation Error for direction array\n");exit(1);
		}
		//Parse the file content into the speed and direction arrays
		int parse_result = WBP_parse(WB_arr_speed, WB_arr_direction);
		if(parse_result == 0){
			printf("File is succesfully parsed\n");
			register int i;
			for(i=0;i<WB_Total_Data_Count;i++){
				*(WB_arr_speed+i) = *(WB_arr_speed+i) / 10.0f;
//				if((i % 24) == 0) printf("\nL(%d) ",(i/24));
//				printf("%g ",*(WB_arr_speed+i));
			}
//			printf("\n%d %d\n",WB_Total_Data_Count,i);
			WB_arr_speed_ext = WB_arr_speed;
			return 0; //Success
		}else{
			printf("\nParsing problem. line = %d\n",parse_result);
		}
	}
	return -1; //Failure
}

/**
 * Function: WB_Finalize
 * ----------------------------------
 * Releases resources used and frees memory allocated by WB.
 * You have to use this routine before exit the application.
*/
void WB_Finalize(){
	if(WB_arr_speed != NULL) free(WB_arr_speed);
	if(WB_arr_direction != NULL) free(WB_arr_direction);
	if(WB_histo_ui != NULL) free(WB_histo_ui);
	if(WB_histo_mi != NULL) free(WB_histo_mi);
	if(WB_histo_pui != NULL) free(WB_histo_pui);
	if(WB_histo_fui != NULL) free(WB_histo_fui);
	WBP_Finalize();
}

void WB_extrapolate_speed(){
	if(WB_height != WB_default_height){
		if(WB_height != WB_height_old){
			if(WB_arr_speed_ext != NULL){
				if(WB_arr_speed != WB_arr_speed_ext){
					free(WB_arr_speed_ext);
					printf("Memory freed 1\n");
				}
			}
			extrapolate_speed(WB_arr_speed, WB_arr_speed_ext, WB_Total_Data_Count,WB_default_height,WB_height,WB_z0);
			printf("Wind speed Extrapolated from %g to %g\n",WB_height_old, WB_height);
			WB_height_old = WB_height;
		}
	}else{
		if(WB_arr_speed != WB_arr_speed_ext){
			if(WB_arr_speed_ext != NULL){
				if(WB_arr_speed != WB_arr_speed_ext){
					free(WB_arr_speed_ext);
					printf("Memory freed 1\n");
				}
			}
			WB_arr_speed_ext = WB_arr_speed;
			printf("Wind speed Extrapolated from %g to %g\n",WB_height_old, WB_height);
			WB_height_old = WB_height;
		}
	}

}

float WB_speed_max(){
	if(WB_speed_max_var>-1) return WB_speed_max_var;
	WB_speed_max_var = get_max_speed(WB_arr_speed_ext,WB_Total_Data_Count);
	return WB_speed_max_var;
}

float WB_speed_min(){
	if(WB_speed_min_var>-1) return WB_speed_min_var;
	WB_speed_min_var = get_min_speed(WB_arr_speed_ext,WB_Total_Data_Count);
	return WB_speed_min_var;
}

float WB_speed_avg(){
	if(WB_speed_avg_var>-1) return WB_speed_avg_var;
	WB_speed_avg_var = speed_avg(WB_arr_speed_ext, WB_Total_Data_Count);
	return WB_speed_avg_var;
}

float WB_speed_median(){
	if(WB_speed_median_var>-1) return WB_speed_median_var;
	WB_speed_median_var = speed_median(WB_arr_speed_ext, WB_Total_Data_Count);
	return WB_speed_median_var;
}

float WB_speed_mode(){
	if(WB_speed_mode_var>-1) return WB_speed_mode_var;
	WB_speed_mode_var = speed_mode(WB_arr_speed_ext, WB_Total_Data_Count);
	return WB_speed_mode_var ;
}

double WB_speed_sd(){
	if(WB_speed_sd_var>-1) return WB_speed_sd_var;
	WB_speed_sd_var = standart_deviation(WB_arr_speed_ext, WB_Total_Data_Count);
	return WB_speed_sd_var;
}

double WB_speed_variance(){
	if(WB_speed_variance_var>-1) return WB_speed_variance_var;
	WB_speed_variance_var = variance(WB_arr_speed_ext, WB_Total_Data_Count);
	return WB_speed_variance_var;
}

double WB_lambert_persistance(){
	if(WB_lambert_persistance_var>-1) return WB_lambert_persistance_var;
	return lambert_persistance_value(WB_arr_speed_ext, WB_arr_direction, WB_Total_Data_Count);
	return WB_lambert_persistance_var;
}

double WB_autocorrelation(){
	if(WB_autocorrelation_var>-99) return WB_autocorrelation_var;
	return autocorelation(WB_arr_speed_ext, WB_Total_Data_Count, 1);
	return WB_autocorrelation_var;
}

double* WB_sd_hourly(){
	double sd_hourly[24];
	standart_deviation_hourly(WB_arr_speed_ext, WB_Total_Data_Count, &WB_height, &sd_hourly[0]);
	register int i;
	for(i=0;i<24;i++) printf("%e\n",sd_hourly[i]);
	return &sd_hourly[0];
}

double* WB_variance_hourly(){
	double var_hourly[24];
	variance_hourly(WB_arr_speed_ext, WB_Total_Data_Count, &WB_height, &var_hourly[0]);
	register int i;
	for(i=0;i<24;i++) printf("%e\n",var_hourly[i]);
	return &var_hourly[0];
}

void WB_variance_monthly(double *p_array_12){
		variance_monthly(&WBP_data[0][0], WBP_Buffer_size, WB_arr_speed, WB_Total_Data_Count, &WB_height, p_array_12);
}

void WB_sd_monthly(double *p_array_12){
	sd_monthly(&WBP_data[0][0], WBP_Buffer_size, WB_arr_speed, WB_Total_Data_Count, &WB_height, p_array_12);
}

void WB_variance_to_sd(double *p_array_variance_12, double *p_array_12,int length){
	register int i;
	/* Calculate the square root of each variance belong to the month*/
	for(i=0;i<length;i++) *(p_array_12+i) = sqrt(p_array_variance_12[i]);
}


void WB_direction_frequencies(int *p_arr_frequency){
	wind_direction_frequencies(WB_arr_direction, WB_Total_Data_Count, p_arr_frequency);
}

void WB_speed_frequencies(double *p_frequency){
	wind_speed_frequencies(WB_arr_direction, WB_arr_speed, WB_Total_Data_Count, p_frequency);
}

void WB_speed_histogram(){
	speed_frequency(WB_arr_speed_ext, WB_Total_Data_Count, WB_histo_ui, WB_histo_mi, WB_histo_pui, WB_histo_fui, &WB_histo_bin_value, &WB_histo_class_count);
}

void WB_SDC(double **p_array_Ti, double **p_array_Wt){
	return SDC(WB_arr_speed, WB_Total_Data_Count, p_array_Ti, p_array_Wt);
}

double WB_energy_pattern_factor(){
	return energy_pattern_factor(WB_arr_speed_ext, WB_Total_Data_Count);
}

double WB_weibull_k_param(){
	if(WB_weibull_k_param_var>-1) return WB_weibull_k_param_var;
	WB_weibull_k_param_var = weibull_k_param(WB_speed_sd(), WB_speed_avg());
	return WB_weibull_k_param_var;
}

void WB_weibull_k_param_monthly(double *k_array_12){
	weibull_k_param_monthly(&WBP_data[0][0], WBP_Buffer_size, WB_arr_speed, WB_Total_Data_Count, &WB_height, k_array_12);
}

void WB_weibull_k_param_by_sd_and_avg(double *p_array_sd, double *p_array_avg, int length, double *p_array_result){
	weibull_k_param_by_sd_and_avg(p_array_sd, p_array_avg, length, p_array_result);
}

double WB_weibull_c_param(){
	if(WB_weibull_c_param_var>-1) return WB_weibull_c_param_var;
	WB_weibull_c_param_var = weibull_c_param(WB_speed_avg(), WB_weibull_k_param());
	return WB_weibull_c_param_var;
}

void WB_weibull_c_param_monthly(double *p_array_avg_speed, double *p_array_k_param, int length, double *p_array_result){
	weibull_c_param_monthly(p_array_avg_speed, p_array_k_param, length, p_array_result);
}

void WB_avg_hourly(double *p_array_24,  double (*operation)(float *, float *) ){
	wind_avg_hourly(WB_arr_speed_ext,WB_Total_Data_Count,p_array_24, &WB_height, operation);
}

void WB_avg_monthly(double *p_array_12, double (*operation)(float *, float *)) {
	Wind_avg_monthly(&WBP_data[0][0], WBP_Buffer_size, WB_arr_speed_ext, WB_Total_Data_Count, p_array_12, &WB_height, operation);
}

void WB_avg_monthly_day_night(double *p_array_12_day, double *p_array_12_night, double (*operation)(float *, float *)){
	Wind_avg_monthly_day_night(&WBP_data[0][0], WBP_Buffer_size, WB_arr_speed_ext, WB_Total_Data_Count, p_array_12_day, p_array_12_night, &WB_height, operation);
}

void WB_avg_seasonal(double  *p_array_4, double (*operation)(float *, float *)){
	Wind_avg_seasonal(&WBP_data[0][0], WBP_Buffer_size, WB_arr_speed_ext, WB_Total_Data_Count, p_array_4, &WB_height, operation);
}

void WB_lambert_persistance_hourly(float *p_array_24) {
	lambert_persistant_value_hourly(WB_arr_speed_ext, WB_arr_direction, WB_Total_Data_Count, p_array_24);
}

void WB_lambert_persistance_monthly(float *p_array_12){
	lambert_persistant_value_monthly(&WBP_data[0][0], WBP_Buffer_size, WB_arr_speed_ext, WB_arr_direction, WB_Total_Data_Count, p_array_12);
}

void WB_lambert_persistance_seasonal(float *p_array_4){
	lambert_persistant_value_seasonal(&WBP_data[0][0], WBP_Buffer_size, WB_arr_speed_ext, WB_arr_direction, WB_Total_Data_Count, p_array_4);
}
//}

//}

